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Abstract 

We introduce a modified SPH approach that is based on discretising the particle 
density instead of the mass density. This approach makes it possible to use SPH 
particles with very different masses to simulate multi-phase flows with large dif- 
ferences in mass density between the phases. We test our formulation with a simple 
advection problem, with sound waves encountering a density discontinuity and 
with shock tubes containing a contact discontinuity between air and Diesel oil. For 
all examined problems where particles have different masses, the new formulation 
yields better results than standard SPH. This is also the case for problems in which 
different spatial resolutions are needed while the mass density does not change. 
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1 Motivation 



SPH is a Lagrangian particle method for solving the equations of hydro- 
dynamics that was invented by Lucy [1] and Gingold and Monaghan [2]. 
Instead of discretising space with a grid, the matter is discretised into so- 
called particles which move with the fluid and do not exchange mass. This 
method is especially suited for compressible flows with irregular bound- 
aries. SPH has been used for many astrophysical problems with great suc- 
cess, and has also been applied to other fields of physics, such as e.g. the 
simulation of liquids [3] and solids [4,5]. 
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Due to its Lagrangian nature, simulating several non-mixing fluids is a 
straightforward extension to SPH. Each particle gets initially marked with 
the phase it belongs to, and these marks do not change with time. 

We are interested in simulating compressible flows with large density dif- 
ferences in contact discontinuities, e.g. due to an interface between a jet of 
liquid near the speed of sound in a surrounding gas. Unfortunately, stan- 
dard SPH as e.g. presented in [6,7] breaks down in this case, as explained 
in this text. Previous approaches to handle large density differences in SPH 
(e.g. [8]) rely on ad-hoc countermeasures. 

For subsonic flows, [9] presents a mesh-free numerical method for gas- 
liquid phase interfaces which is based on the MPS method [10]. It is re- 
stricted to incompressible flows, whereas we are interested in compressible 
flows. We are also interested in similar spatial resolutions in both phases of 
the flow. This requirement is different from e.g. astrophysical collapse sim- 
ulations, where one usually wants a higher resolution in the denser regions. 

Another point of interest for us is having differing spatial resolutions, i.e. 
particles with different masses, in regions where the mass density does not 
change. This has been used on several occasions (e.g. in [11]). However, this 
can lead to rather large errors with standard SPH. 

We will in the following present a modification of SPH that interprets cer- 
tain numerical quantities in a different manner [12], leading to a stable, ro- 
bust, and accurate evolution even in these cases. 



2 Describing multi-phase flows 

One commonly used way of introducing the SPH discretisation (see also [7]) 
starts out by considering an arbitrary field /(x). This field is first smoothed 
by folding it with a kernel W(x), which leads to the smoothed field (/(x)) 



where the kernel W(x) must be normalised according to / d 3 x W(x) = 1. 
One usually chooses kernels that have approximately the shape of a Gaus- 
sian, and that have compact support for reasons of efficiency. The size of 
the domain of support is called the smoothing length and is usually denoted 
with the letter h. 

In the next step, the smoothing integral is discretised at N particle positions 
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Xj, which can in principle be chosen freely, but should of course be "reason- 
ably" distributed. This leads to the SPH approximation f(x) 



N 



/« == I Vj fj W(x - Xj) 



(2) 



of the field. The volumes Vy are the discrete counterparts of the volume 
element d 3 x' in the integral above, and have to be chosen so as to be consis- 
tent with the spatial distribution of the particles. It is customary to assign a 
certain mass rrij to each particle, and then replace the volumes through 



where p ; is the discretised mass density assigned to the particle. This is 
motivated by the fact that the particles do not exchange mass, which leads 
to the time evolution equation 



making m z a natural choice for one of the primary variables. In order to 
make SPH Lagrangian, the particles have to move with the flow, leading to 



as the time evolution equation for the particle positions Xj. Here Vj is the 
discretised fluid velocity field. 

Assuming that the desired spatial resolution is about the same in all fluid 
phases, one would choose similar particle spacings there, leading to similar 
particle volumes Vj. If the mass densities in the different fluids are about the 
same, no further problems arise, and SPH as usual can be used to describe 
them. However, if the difference in mass density is large (say, about one 
order of magnitude or more), then the particle masses will differ by the 
same factor, leading to problems at the phase interfaces. These problems 
are mostly caused by inaccuracies in the mass density. These inaccuracies 
are substantial for largely different particle masses. Below we describe how 
they come about, and how they can be avoided. 
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2.1 Standard SPH 



There exist in principle two different methods for obtaining the discretised 
mass density p\ in the "standard" SPH formalism. Both start out by con- 
sidering the approximate mass density p at the particle position x/. This 
quantity is obtained from eqn. (2) by choosing /(x) = p(xf) and applying 
eqn. (3), leading to 

p(x l )=£m ; W u - (6) 

i 

where the abbreviation W;y := W(x ; — xy) has been introduced. 

The first and conceptually simpler method to obtain p\ from this is by set- 
ting pi := p(xf), leading to 

Pl = £mjW tj . (7) 

;' 

This method is often used for astrophysical problems when there are free 
boundaries, i.e. when the matter distribution extends into vacuum. How- 
ever, it is not suited for a phase interface with a large density discontinuity. 
The smoothing inherent in eqn. (7) will smooth out the density jump over 
a region of the size of the smoothing length h in either direction of the dis- 
continuity. Particles in this region will then "see" a density that is much 
less or much larger from the real density in their phase. When this density 
is used to calculate the pressure through the two phases' equations of state, 
the pressure will be very wrong (as shown in figure 1 for two ideal gases), 
and it is basically impossible to set up a stable interface in equilibrium. This 
problem becomes even more severe when one of the fluids has a stiff equa- 
tion of state, as is the case e.g. in liquids or solids, because then density 
inaccuracies will lead to even larger errors in the pressure. 

The second method for obtaining the mass density pi is by integrating pi in 
time via the time derivative of eqn. (7), leading to 



^ P! = £m ; (v ; -v ; ).VW !7 (8) 

where eqn. (5) has been used, and the abbreviation VW/y := ( VW) (x, — xy) 
has been introduced. This method has the advantage that the initial data 
for Pi can be chosen freely, so that density discontinuities can be modelled. 
This can be used to simulate surfaces of liquids and solids. 
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Figure 1. Physical density and pressure and the corresponding standard SPH quan- 
tities at a phase interface with equal pressure and a density ratio of 100 : 1. Near 
the interface, the approximation errors reach a factor of about 2 in the dense re- 
gion and about 50 in the thin region. Essentially the same picture results when 
there is no change in mass density but when instead the particle mass and spacing 
are changed by a factor of 100 to achieve a higher spatial resolution (not shown). 
(These approximation errors do not exist when a different SPH formulation is used; 
see section 2.2 below). 

The problem with this method is similar to the problem encountered when 
calculating the density directly from the particle distribution by eqn. (7). 
The particles on each side of a phase interface "see", via the term ray in eqn. 
(8), very different particle masses on the other side of the interface. The 
values of d/dt Pi then contain large inaccuracies, leading to instabilities at 
phase interfaces. 



2.2 Modified SPH 



However, eqn. (7) is not cast in stone. An ansatz equivalent to but different 
from the one leading to this equation is not to smooth the mass density p,-, 
but rather the particle density n\ = 1/Vj [12]. This is easily motivated by 
the fact that the mass density can be discontinuous over a phase interface, 
while the particle density is not, according to our assumption of similar 
spatial resolutions on both sides. Smoothing the particle density n\ := n(xf) 
via eqn. (2) leads to 



(9) 



and by taking its time derivative, the equation 



£«i = E(v,-v y )-VW f ,- 



(10) 
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is obtained after using eqn. (5). As it is customary in SPH to use Pi instead 
of tii, we apply eqns. (3) and (4) and arrive at 

j t Pr = m 1 ^(v l -Vj)-VW l j . (11) 

This new formulation of the equation of continuity is the key element of 
our SPH approach. It should be noted that this equation is identical to eqn. 
(8) when all particle masses nij are the same, which is the case for many 
single-phase SPH simulations. 

For the simulations presented in this text, we discretise the Euler and the 
internal energy equations in established ways: 

d 1 Ttl ' 

S v, = --E^(P, + P,) VW, (12) 
d 11 Ttl ' 

5«=2 5 E^fa + *)ta-'y)-™ (/ d3) 

where e\ is the specific internal energy. The symmetrisations (py + p,) are 
e.g. explained in [7]. 

3 Tests 

In the following we test the new SPH formulation and compare it to an- 
alytic solutions as well as simulations performed using standard SPH as 
described in [7]. That is, the only difference between these two formula- 
tions is that we use eqn. (11) instead of eqn. (8). This also means that other 
parts of an SPH code such as time integration or artificial viscosity are not 
affected. As test cases, we use an advection problem, a sound wave encoun- 
tering a discontinuous change in the sound speed, and a shock tube with a 
Diesel-air interface. 



3.1 Advection equation 

We compare the standard and the new SPH formulation by simulating a 
one-dimentional advection equation. That is, we solve the equation of con- 
tinuity for the density p while prescribing the velocity field v. The velocity 
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Figure 2. Comparison of simulations with the standard SPH and the new SPH 
formulation at five different times. The dotted lines show the analytic solution, the 
solid lines show the particle values. Standard SPH is less accurate at the narrower 
peaks, which correspond to later times. 

profile (which is constant in time) and the initial density profile are given 
by 



*M = ITj? (14) 

Po (x)=Ax 2 exp{-(^p) 2 ) (15) 



with A — 1.5, xq = 1, W — 0.4, and q = 0.2. We initially place the parti- 
cles with equidistant spacings with a density of n = 10 particles per unit 
length and use a smoothing length of h = 0.25. The particle masses are cho- 
sen according to the initial density profile at the initial particle positions, 
i.e. they differ. Advection problems are particularly well suited test prob- 
lems for Lagrangian methods, so we expect a high accuracy from this low 
resolution. 

The results of simulating this equation with both the standard and the new 
SPH formulation are presented in figure 2, which shows the density p at five 
different times. Both formulations track the analytic solution very nicely in 
spite of the coarse resolution. However, at later times, the standard SPH 
formulation underestimates the density near the peaks, while the new for- 
mulation stays much closer to the analytic solution. 



3.2 Sound wave 



The sound wave test case consists of two regions containing the same ideal 
gas, but with different densities and in pressure equilibrium. The density 
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Figure 3. A sound wave crossing an interface between two ideal gases with a den- 
sity ratio of 10 : 1. The interface is at x = 0. Shown is the pressure at the times 
t = —0.5, t = 0, and t = +0.5, where dotted lines mark the analytic solution. The 
graph in the lower right hand corner shows the result at t = +0.5 of a simulation 
with standard SPH for comparison. 



discontinuity is located at x = with a density ratio of 10 : 1. These condi- 
tions also lead to different temperatures, and to sound speeds with a ratio 
of 1 : \/W. Figure 3 shows an initially Gaussian-shaped sound wave at 
different times. 



At t = —0.5 the initial wave travels to the left. At f = the wave has 
reached the interface where it is partially transmitted and partially reflected. 
At t = +0.5 the wave consists of two packets, travelling in different direc- 
tions with different speeds. The simulation was performed with n = 200 
particles per unit length and a smoothing length of h = 0.1. The analytic 
solution is shown as dotted line underneath the simulation result. The SPH 
simulation with the new equation of continuity (11) tracks the analytic so- 
lution quite well. On the other hand, standard SPH using eqn. (8) performs 
rather poorly in this case, as can be seen in the graph in the lower right 
hand corner: the transmission and reflection coefficients are wrong, and the 
pressure develops spikes at the interface. We assume that the reason for this 
is just the one demonstrated in figure 1. 
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Figure 4. A shock wave emanating from an interface with air to the left and Diesel 
to the right. Shown are pressure and density at t = 5 x 10~ 4 s, where dotted lines 
mark the analytic solution. The discontinuity is initially at x = 0. The shock front 
is spread out over 8 smoothing lengths by artificial viscosity. The pressure differ- 
ence across the rarefaction wave is rather small. The shock and rarefaction waves 
are not visible in the density graph because of the large scale differences. The con- 
tact discontinuity is stable and well preserved in spite of the coarse resolution; the 
pressure oscillations do not grow with time. 



Table 1 

Initial data for the shock tube test case 



Quantity 


air 


Diesel oil 


P [kg/m 3 ] 


8.81 


772.546 


p [MPa] 


10 


5 


T [K] 


3931.5 


393.15 



3.3 Shock tube 



A further test case for our formulation is a shock tube containing a Diesel- 
air interface, shown in figure 4. The initial discontinuity is at x = 0, with 
air to the left and (liquid) Diesel oil to the right. The initial pressure ratio 
is 2 : 1, the density ratio about 1 : 80. Here we have artificially decreased 
the air density and increased the air temperature by a factor of 10 to create 
a more difficult test case. Table 1 lists the exact initial data for this test case. 
The shock wave in the Diesel oil travels to the right, the rarefaction wave in 
the air to the left. Because the Diesel oil is nearly incompressible, the final 
pressure is close to the initial pressure in the air phase. The equation of state 
for the Diesel oil was kindly provided to us by the Robert Bosch GmbH. 

The simulation was performed with n = 100 particles per metre and a 
smoothing length of h = 5 x 10" 2 m. We use the artificial viscosity pre- 
sented in [7] with a viscosity coefficient of a = 0.5, because some artificial 
viscosity is necessary to produce entropy in the shock front. This spreads 
out the shock front over about 8 smoothing lengths, which is acceptable for 
our purposes. 
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Table 2 

Comparison of the analytic solution and the simulation result for the shock re- 
lations in the Diesel phase. v s is the shock front speed, v D the post-shock Diesel 
speed. The sound speed in the pre-shock Diesel phase is 1059.6 m/s. 



Quantity 


analytic 


SPH 


v s [m/s] 


1077 


1056 ± 20 


v D [m/s] 


5.93 


5.94 ± 0.02 


Ap D [kg/m 3 ] 


4.27 


4.28 ± 0.02 


Ap D [MPa] 


4.93 


4.94 ± 0.02 


AT D [K] 


0.860 


0.855 ± 0.005 



The initial pressure discontinuity remains visible as spikes at the contact 
discontinuity. These spikes are caused by the numerical initial data, which 
have a discontinuity and hence contain high frequency modes that are not 
resolved in the simulation. According to eqn. (1), the initial data should be 
smoothed before the SPH formalism is applied. We skip this step because 
we want to show that these high frequency modes do not harm the simula- 
tions. They remain present, but are not amplified. The formulation is stable. 
Table 2 compares several important quantities of the simulation result to the 
analytic solution, showing very good agreement in the shock relations. 

We did not manage to perform this simulation with standard SPH. As il- 
lustrated in figure 1, the error in the density near the contact discontinuity 
leads to an error in the pressure. The pressure error can be orders of magni- 
tude larger than the pressure itself for fluids which have a stiff equation of 
state. This problem does not exist in the new formulation. 



4 Conclusion 

We describe a modification to the standard SPH formalism that smoothes 
the particle density instead of the mass density. As tested by simulating an 
advection equation, sound waves, and shock waves, the new formulation 
either yields more accurate results than standard SPH, or the equivalent 
simulation with standard SPH is not stable. We conclude that this modified 
SPH formulation is an effective method for simulating multi-phase flows, 
and conjecture that this modification is beneficial for all simulations where 
particles have different masses. 

The authors would like to thank Prof. Hanns Ruder for encouraging this 
work. We would also like to thank the Robert Bosch GmbH for providing 
an equation of state for Diesel oil. Financial support was provided by the 
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